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Abstract 

A revised and greatly improved version of the 3D continuum radiative transfer 
code MC3D is presented. It is based on the Monte-Carlo method and solves the 
radiative transfer problem self-consistently. It is designed for the simulation of dust 
temperatures in arbitrary geometric configurations and the resulting observables: 
spectral energy distributions, wavelength-dependent images, and polarization maps. 
The main objective is the investigation of "dust-dominated" astrophysical systems 
such as young stellar objects surrounded by an optically thick circumstellar disk 
and an optically thin(ner) envelope, debris disks around more evolved stars, asymp- 
totic giant branch stars, the dust component of the interstellar medium, and active 
galactic nuclei. 

PACS: 95.30.Jx - Radiative transfer 

Key words: Continuum radiative transfer; Monte-Carlo method; Numerical 
simulation; Polarization; Circumstellar shells 



1 Introduction 



The 3D continuum radiative transfer (RT) code, MC3D, combines the most re- 
cent Monte Carlo (MC) radiative transfer concepts for both the self-consistent 
RT, i.e., the estimation of spatial dust temperature distributions, and pure 
scattering applications, taking into account the polarization state of the radi- 
ation field. It has been tested intensively and compared with grid-based and 
MC RT codes (see, e.g., Pascucci et al. (|l]). Wolf (|)). 

In addition to the previous version of MC3D (Wolf & Henning (^; see also 
Wolf et al. (|^)), MC3D (V2) contains several improvements which allow a much 
more efficient solution of the radiative transfer problem in three-dimensional, 
arbitrary dust configurations. The main features are: 
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• A new concept for estimation of dust temperature distributions which was 
first described by Bjorkman & Wood (|^) in the context of the solution of the 
radiative transfer problem in circumstellar disks. It is based on the imme- 
diate correction of the dust grain temperature and the reemission radiation 
field after absorption of a photon packet. Beyond a faster estimation of the 
dust equilibrium temperature without iterations, it allows the application 
of the code for optical depths far above the previous r ^ 10^ limit. 

• A new concept for the test photon transfer, based on the work of Lucy 
(previously developed for one-dimensional model geometries), which takes 
into account the absorption not only at the points of interaction of the test 
photons with the dust grains but also in between. Therefore, the tempera- 
ture correction can be performed in all volume elements on the test photon 
path and consequently, there is no more lower optical depth limit for the 
application of MC3D. 

• A new concept for the simulation of scattering in optically thin dust / electron 
configurations which has been firstly described by Cashwell & Everett (|^) 
(see also Witt (|^) and Menard (P) for former applications in one- and two- 
dimensional model geometries) which makes use of the enforced scattering 
method. 

• The simulation of the RT in dust grain mixtures consisting of grains with 
different chemical composition and size. 

• A raytracer for the simulation of images and spectral energy distributions 
(SEDs) which may be used for the simulation of the dust reemission in 
the midinfrared-millimeter wavelength range, where scattering is of minor 
importance. 

• An accelerated solution of the temperature estimation in flat disks with or 
without vertical density structure. 

These features have been adapted to and tested for several ID, 2D, and 3D 
model geometries which are shown in Fig. 1. 

Previous applications of MC3D cover feasibility studies of extrasolar planet 
detections (Wolf et al. (0)), the RT in the clumpy circumstellar environment 
of young stellar objects (Wolf et al. (pJ])), polarization studies of T Tauri stars 
(Wolf et al. QT^) ), AGN polarization models (Wolf & Henning (|13D), a solution 
for the multiple scattering of polarized radiation by non-spherical grains (Wolf 
et al. (0)), and the inverse RT based on the MC method (Wolf ([T|)). 

The executables of MC3D (V2) can be downloaded for several model geome- 
tries and computer platforms from 

[http : //www ■ mpia-hd . mpg . de/FRINGE/SDFTWARE/inc3d7 

(current US mirror page: 

[http : //spider . ipac . caltech. edu/staf f /swolf /nic3d/[ ). 
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Fig. 1. Model geometries: Examples for the subdivision of the model space into 
energy storage cells (see Sect. 3.3 for further explanations). ID: radial symme- 
try; 2D (a): radial symmetry inside a disk; 2D(b): fully two-dimensional model 
for configurations with a radial and vertical dependence of the density distribu- 
tion; 3D(a-c): three-dimensional models considered in cartesian, cylindrical, and 
spherical coordinates. 



Furthermore, a detailed description, usage instructions, and IDLQ routines 
for a subsequent data analysis of the code are provided there. Very frequently 
used model geometries, such as disks with a radial and vertical variable den- 
sity distribution, which is often used for the analysis of SEDs, images and 
polarization maps of circumstellar disks (see, e.g.. Wood et al. (0), Cotera 
et al. ([T7|) , Fischer et al. ([I8|)), are available. However, the author encourages 
the reader to ask for further model geometries. 



In Sect. 2 the general scheme for the solution of the RT problem on the basis of 
the MC method as it is applied in MC3D is briefly outlined. In Sect. 3 the two 
different concepts for the estimation of the spatial dust temperature distri- 
bution (iterative radiative transfer and immediate reemission) are discussed. 
In the following Sect. 4, the determination of the observables is described, 
and Sect. 5 introduces further concepts which are used in MC3D in order to 
increase the efficiency of the RT algorithms. 



Interactive Data Language 
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2 General remarks about the applied solution of the radiative 
transfer problem 

Each model to be considered with MC3D consists of two independently defined 
components: 

(1) The spatial density distribution of the scattering and absorbing medium, 
and 

(2) The radiative source(s). 

Both must be defined inside a convex model space in order to exclude the case 
of radiation leaving the model space at one point and entering it at another. 



2.1 Scattering and absorbing medium 

MC3D handles spherically symmetric scatterers/absorbers, such as spherical 
dust grains and electrons (in the Thomson scattering regime). The required 
optical properties, i.e., the absorption and extinction efficiency (Qabs, Qext), 
the special Mueller matrix S (see Sect. 2.3, App. C), and quantities to be 
derived from these - the albedo A, the extinction and scattering cross section 
(Cext, Qsca), and the scattering function - can be calculated with an embed- 
ded Mie scattering routine (using the Mie scattering algorithm published by 
Bohren & Huffman (0)) on the basis of the real and imaginary refractive 
index (n,fc) and the particle radius a. Since Thomson scattering can be - 
formally - considered as a special form of the Mie scattering formalism (no 
absorption, no wavelength dependence), in the following only dust grains are 
considered as the scattering and absorbing medium, covering also the special 
case of electrons. 

MC3D allows simulation of the RT in arbitrary density distributions Pj{r), 
where r represents the spatial coordinate and the index j refers to the consid- 
ered dust species (characterized by the grain size and chemical composition). 
If MC3D is used to simulate the scattering by dust grains only (for instance for 
the simulation of images, polarization maps or SEDs in the visual-midinfrared 
wavelength range), the mean optical properties of dust grain mixtures con- 
sisting of 3 chemically different components and a size distribution following 
a power-law n{a) oc a~" {n{a) is the number of dust grains with the radius 
a; a=const.) can be estimated on the basis of the formalism being described 
in Appendix A. As found by Wolf (pOD , this concept of mean dust grains also 
represents a reasonable approximation for the estimation of temperature dis- 
tributions in optically thick density distributions and the simulation of SEDs 
in general. 
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The density distribution may be defined either analytically or on a predefined 
grid in order to allow the implementation of density distributions resulting 
from hydrodynamical simulations (for examples, see Wolf et al. Sect. 4; 
Wolf et al. (|T0D). 



2.2 Radiation sources 



In principle, radiation sources, that are arbitrary with respect to their num- 
ber and spatial configuration, spatial extent, intrinsic SED, and radiation 



characteristic can be considered (see Wolf et al. (P, (|T0D for examples). In 
MC3D the RT is simulated at certain wavelengths within the wavelength 
range [Amin, Amax]- For this reason, the monochromatic luminosity[3] Lx of the 
radiation source(s) is partitioned into nphoton so-called "weighted photons" 
(Sobolev (0)) each of which is characterized by its wavelength A and Stokes 
vector / = (/, Q, U, V)'^. The polarization state of a photon is determined by 
the Stokes vector components as follows: 



(a) Linear polarization degree : Pi = V*^*^^"^^ I 

(b) Orientation of the linear polarization : 7 = | arctan (^^^ 

(c) Circular polarization degree : Pc = j 



The two main radiation sources in MC3D RT simulations shall be introduced 
in brief: 

(1) Pointlike Star: Here the emission direction of photons, described in spher- 
ical coordinates (6'e, 4>e) is given by 

COS^E = -l + 2Zi, 0E = 271^2, (2) 

(2) Extended Star (radius Rs, radiation characteristic at each point on the 
stellar surface: / oc cos(6'e)): 

sin 6'e = \/^, 0E = 2 vr ^4 , (3) 

whereby 9'-^ and 0^ are related to the z' axis parallel to the radius vector 
{Rg, 6'e, 0e) which can be determined using Eq. 2 (see Fig. 2 for explana- 
tion). Here and in the following, Zi {i = 1,2, . . .) are random numbers uni- 
formly distributed in the interval [0,1]. For their determination a random 
number generator, which minimizes the sequential correlations between ran- 
dom numbers by the combination of different random number generators, is 

^ Between the (bolometric) luminosity L and the monochromatic luminosity Lx the 
following relation exists: L = /q°° Lx d\. 
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Fig. 2. Emission angles, [a] Pointlike star, [b] Extended source. 

used (Knuth (p2D). The Stokes vector of a newly emitted, unpolarized photon 
is defined as Jq = (/, Q, U, V)^ = (1, 0, 0, 0)^. 

Assuming the dust grains to be spherical, the direction of photon reemission 
can be determined as for point-like stars shown in Fig. 2 [A] (Eq. 2). Since the 
point of reemission inside certain volume elements of the model space (energy 
storage cells, ESCs - see Sect. 3.3) is assumed to be constant, the location 
of the reemitting dust grains has to be chosen randomly therein, assuming a 
constant density distribution within the cell. The determination of the point 
of reemission within an ESC is described for the different model geometries 
shown in Fig. 1 in App. D. 



2.3 Photon transfer and scattering 



The mean free path length / between the point of emission and the point of 
the first photon-dust interaction (and - subsequently - between two points of 
interaction) is given by (see App. B) 



Toxt,i = -ln(l 



(4) 



^cnd 

"^ext,! = X! 

i=l 



no 

^Pj(ni) -Cextj 
i=i 



1=1 



(5) 



Here, ud is the number of different dust grain species and f\. is the spatial 
coordinate corresponding to the ith integration point along the path length /. 

The step width A/j has to be small enough to satisfy the linear approximation 
of the extinction cross section and density distribution along the path of the 
photon (Eq. 5). The step width A/j is not fixed but has an initial (maximum) 
value which is chosen to be a fraction of the extent of the current energy storage 
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cell (ESC, see Sect. 3.3 for an introduction to the ESC concept) the photon is 
moving through - or, if any of the neighboring cells is smaller, then of that. 
This procedure ensures that no ESC is skipped on the photon's path. Because 
the subdivision into ESCs is adapted to the model geometry and distribution 
of the local optical depth, this coupling drastically increases the efficiency of 
the photon transfer since a small step width in high density regions as well as 
large step widths in low density (gradient) regions is provided simultaneously. 
In order to minimize the difference between the exact value rext,i and the 
numerically achieved value fext.b is decreasing as soon as the test photon 
reaches the vicinity of the next interaction point. 

This photon transfer concept was found to limit the applicability of previous 
version of MC3D to optical depths r > 10~^ (Wolf & Henning (|^)) since the 
probability for photons to interact with the (circumstellar) matter is negligibly 
small for even lower optical depths (it decreases exponentially as described by 
Bouguer-Lambert's law). In order to be able to consider lower optical depths, 
for instance the simulation of the scattered light arising from an optically 
thin envelope or the electron environment around hot stars, the concept of 
enforced scattering (Cashwell & Everett (|^)) has been implemented. The idea 
is to force each photon to be scattered at least once between the point of 
emission and the boundary of the model space, provided there is dust on 
its path (optical depth Texto)- A fraction of the photon (e~'^™*o) leaves the 
model space without interaction while the remaining part (1 — e~'^°'"*o) will be 
scattered. Following the same strategy as for the determination of the mean 
path length (see App. B), we derive: 

t (rext)drext - . drext " 1 - e---o 



Z = F(rext) = / /(Text) dr'xt = 





Text = -ln(l-Z[l-e-"^-o]). (8) 



As Fig. 3 illustrates, the enforced scattering concept is of tremendeous impor- 
tance for the simulation of scattered light of, e.g., the optically thin envelope 
of young stellar objects, debris disks aroung more evolved stellar systems, or 
the zodiacal light in extra-solar planetary systems. 

Depending on the applied concept for the temperature estimation (see Sect. 3), 
scattering or absorption or both processes occur at a point of interaction. The 
modification of the Stokes vector due to the ith scattering is described by a 
special Miiller matrix SjiO^cj)) of the dust grain species number j (see, e.g.. 
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Fig. 3. Illustration of the high efficiency of the enforced scattering concept in opti- 
cally very thin configurations, (a) Without and (b) with enforced scattering. Spher- 
ical shell being illuminated by an embedded star. Density profile: p{r) oc r~^'^. 
Optical depth at the considered wavelength (as seen from the star): 10~^. Number 
of photons: 10^. Shown is the (Intensity) without the central star. The required 
CPU time amounts in both cases to about 7min (using a Intel Celeron, 1.2GIIz 
processor). 

Bickel & Bailey (|3|), Bohren & Huffman (|19|)): 

i,^s,{e,<p)fi_, . (9) 

Here, 6 and are the scattering angles (see App. C). The estimation of the 
scattering direction, based on the scattering function provided by Mie calcula- 
tions, and the consideration of multiple scattering events, taking into account 
(the change of) the polarization state of the weighted photon, are based on 
the concepts described by Fischer et al. (p^, App. A3 and A 4.2. The absorp- 
tion and reemission processes are described in the context of the temperature 
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estimation in Sect. 3. 



3 Estimation of the spatial temperature distribution 



Two different concepts for the solution of the self-consistent RT problem, i.e., 
the determination of the spatial dust temperature distribution on the basis of 
the amount of absorbed energy of the stellar and the dust reemission radiation 
field, are embedded in MC3D. While the first is based on an iterative procedure 
(Sect. 3.1), in the second the temperature distribution is changed ("updated") 
simultaneously with each absorption act (Sect. 3.2). 



3.1 Iterative Radiative Transfer 



Assuming a star with an effective temperature T^s and radius i?*, its monochro- 
matic luminosity (assuming blackbody radiation) can be written as 



A,Ster: 



-47ri?2.7r5A(reff), 



(10) 



where B\{T) is the Planck function. Due to absorption of stellar and reemit- 
ted radiation the dust is heated. Absorption occurs simultaneously to the 
scattering whereby the probability for a dust grain of the species j to be the 
interaction partner of the photon is given by 



ni(j) 



Pi(r) -Ccxtj 



i'=i 



The difference A/ of the Stokes parameter / before and after the absorption 
process is determined by the wavelength-dependent albedo of the dust grains 



A,/ = /,_i-(l-yl;,.). 



:i2) 



Consequently, the amount of monochromatic luminosity represented by the 
particular photon is decreased by 

A,L;, = • A,7 . (13) 

^Photon!, A j 

The monochromatic luminosity of a single spherical dust grain with the radius 
a, the temperature Tg, and the wavelength-dependent absorption efficiency 



9 



Q^^^ can be written as 

LX^=4na'-Qf^^{a)-nB,{T,). (14) 



Assuming a static dust configuration, the energy being reemitted by all grains 
of the dust species j within the volume element V during the time interval 
At has to be equal to the energy being absorbed by the grains (local energy 
conservation) : 

oo oo 

At ■ Ny. J LI. dX^At-J Lf^^ dA , (15) 





where Ny. is the number of dust grains (species j) in the volume V. The ab- 
sorbed monchromatic luminosity Lf°^^ within the volume V is the accumulated 
amount of absorbed photon "fractions" : 



^fv = E E (a.^aJ, , (16) 

k=l i=0 



where UahsU, k) is the number of absorption acts of the kth photon by the 
dust species j in the volume V. The combination of Eq. 14 and 15 results in 
an equation which allows determination of the mean dust grain temperature 



Tg. in the volume V: 



/o°° Lf' dA 7 

jQt;{a)-B,{f,^)dX. (17) 



Since the temperature Tg. cannot be separated from Eq. 17, in MC3D it is de- 
termined from a table of precalculated values of the right-hand integral in the 
range of K to the dust sublimation temperature. A temperature step width 
of 0.1 K was found to be sufficiently small for most astrophysical applications. 

To consider both, the stellar heating and the subsequent dust reemission ra- 
diation, the following iterative procedure is applied: 

(1) Heating by the star. 

(2) Reemission by the dust grains, based on the energy being absorbed from 
the stellar radiation field. The determination of the reemission coordi- 
nates in the ESCs depends on the configuration and is described in App. D 
in detail. 

(3) Reemission of the dust grains, based on the sum of the energy being 
absorbed from the stellar and the dust reemission radiation field. 
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The third step has to be repeated as long as the difference between the input 
energy (given by the stellar radiation field) and the output energy (given by 
the sum of the attenuated stellar radiation field and the reemitted radiation 
outside the model space) is larger than a user-specific threshold (e.g., 0.1% of 
the input energy). 



3.2 Immediate Reemission 



The second concept embedded in MC3D for the estimation of the spatial dust 
temperature distribution has been developed by Bjorkman & Wood (H). In 
contrast to the iterative concept (Sect. 3.1), either scattering or absorption 
occurs at a point of interaction of a photon with a dust grain. The probability 
for a photon to undergo the one or the other interaction process is given by 



n..(j,ri 



where 'x' stands for either absorption or scattering. When a photon packet 
is absorbed, the new dust temperature is calculated immediately based on 
Eq. 17. To correct the ESC temperature, the wavelength of the immediately 
reemitted photon is chosen so that it corrects the temperature of the spectrum 
previously emitted by this particular ESC, i.e., the probability distribution of 
the reemission wavelength is given by 

/(A) (X Qabs [5A(Tg.(^)) - B^{f^.{i - 1))], (19) 

where T^.^{i — 1) represents the temperature of the grain species j before and 
Tg. (z) after the zth absorption process in the considered ESC. The value of 
the monochromatic luminosity of the reemitted photon results from the cor- 
responding energy conservation equation: 

L^,d\^ = Lx,dX2, (20) 



where dXi and d\2 are the wavelengths of the absorbed/reemitted photon. As 
the simulation runs, the temperature distribution (and the SED given by the 
photons leaving the model space - see Sect. 4) relaxes to its equilibrium value 
without iterations. 

To reduce the number of photons being required for the temperature estima- 
tion in optically thin configurations, the concept of Lucy (|^) is apphed. It 
considers the absorption of the electromagnetic radiation field not only at the 
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end points of the photon path (points of interaction) but also in between. The 
absorption due to aU dust grain species has to be taken into account: 



no ^"7 



point2 



''"abs — 



Cabsj • dr. 



(21) 



^ point 1 



The example of a circumstellar disk with an optically thin envelope shown in 
Fig. 4 illustrates the great potential of this concept. 



3.3 ESC concept 

To achieve a spatially resolved dust temperature distribution, the model space 
is subdivided into so-called energy storage cells (ESC) in which the absorbed 
energy of the radiation field is stored (accumulated) and in which a constant 
dust temperature is assumed. To represent the spatial temperature structure 
according to the model geometry, defined by both the spatial distribution of 
radiation sources and the dust density, the model space has to be subdivided 
taking this into account. Furthermore, the spatial distribution of the optical 
depth has to be taken into account in order to resolve the high temperature 
gradient in optically thick regions. Finally, the size of the ESCs also has to 
decrease towards the radiation source even in case of optically thin configura- 
tions to trace the steep temperature gradient resulting from the geometrical 
dilution of the radiation field. Fortunately, the last two conditions are fulfilled 
simultaneously for most astrophysical applications since the density gradient 
increases towards the radiation source(s). 

While the ESC distributions ID, 2D(a), 2D(b), 3D(b), and 3D(c), shown in 
Fig. 1, arc fully determined by user specific parameters (such as the number of 
ESCs and the refinement factors for the model space towards the central object 
and/or the disk midplane), model 3D (a) is prepared to provide an adaptive 
ESC generation according to the constraints described above. The model space 
described by model 3D(a) is a cube with the side length Ic- It is subdivided 
according to the distribution of the optical depth. Here, the quantity 



is defined, where Ny{j) is the number of dust grains (species j) in the volume 
element Y {V < l^). At the first step, pv is estimated inside the whole model 
space (pv — Po)- In the second step, the model space is subdivided into 8 



PV = E^ext(j)-^v(j) 



(22) 
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cubes (side length Zc/2). In each "sub-cube" the relation 

^<S<1 (23) 
Po 

is checked, whereby pv(^) is the value of py in the cube and H is a user- 
specific quantity which controls the subdivision depth (number of subdivision 
levels). In the subsequent steps all sub-cubes for which the relation (23) is not 
fulfilled, the subdivision is continued. Finally, the model space is subdivided 

into cubes with the side length " 1^ {ic = 1,2, . . .). To prevent this algo- 
rithm ending up with too large a number of cubes (large RAM requirement), 
resulting from a too small value given for S, the maximum number of cubes 
is a second user-specific parameter (constraint) of this algorithm. If the num- 
ber of cubes exceeds this value, S is increased, the subdivision process starts 
again, and so on. 



4 Observables 



Polarization maps, images, and SEDs are the obervables which can be directly 
simulated with MC3D. For this reason the RT transfer is simulated indepen- 
dently for both the star and the dust as the sources of emission. The RT is 
performed as described in Sect. 2.3 and Sect. 3.1, Eq. 10-14 but without it- 
erations. Although the immediate reemission concept provides the SED and 
images as a direct result of the RT too (see Sect. 3.2), it is not discussed here 
since the polarization state of the radiation is not considered and a decou- 
pling of the short-wavelength stellar radiation and the long-wavclcngth dust 
reemission (and the simultaneous decoupling of the radiation at separate wave- 
lengths) provides more freedom for the selection of the considered wavelength 
(range) for which the observables shall be simulated. See also Sect. 5 for a 
description of the implemented raytracer which can be used to derive the dust 
reemission SED and images. 

In order to derive the observables, a photon will be projected onto an observing 
plane oriented perpendicular to the path of the photon as soon as its next 
scattering position would be located outside the model space (see Fig. 5 [A]). 
Consequently - in approximation of real observations - only those photons 
will be projected onto the same plane which leave the model space in parallel 
directions. The number of observing planes is - with regard to the possible 
last scattering direction - potentially infinite. 

The next step is to combine all those observing maps to a single map for which 
the angular distance of their normal vector to a given direction is smaller 
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than the angle uj (see Fig. 5[B]). The sum of the monochromatic luminosity 
of all photons collected on this map therefore represents the monochromatic 
luminosity L\ q^ {L\ n < L'^J of the object seen under the solid angle = 
27r(l - cosa;[rad])- 

Let (Ja be the radiating area of an object seen in projection on the observing 
plane. Then, the corresponding monochromatic intensity Ix of the object is 

h = ^ . (24) 

The flux density S\ at the point of observation is determined by the distance 
d between the object and the observer as follows: 

Sx = ^-h. (25) 

In order to derive spatially resolved images (such as shown in Fig. 4), the ob- 
serving planes are subdivided in pixels. Regarding the last scattering position 
of a photon, its (projected) Stokes parameters, wavelength and monochromatic 
luminosity are stored in the corresponding pixel. The required transformation 
of the Stokes vector was outlined by Fischer (|25|). 



5 Diverse further embedded concepts 

(1) Selective Emission (for self-consistent RT) 

(a) Wavelength range selection: 

The CPU time can be drastically reduced by considering only those 
wavelength ranges in which the particular radiative source (the star 
or the dust in a certain ESC) is emitting a non-negligible amount of 
energy (see Wolf & Henning 2000) for details. 

(b) Geometrical selection: 

For certain astrophysical models, the efficiency can be improved by 
excluding photons which would be emitted into dust free regions. 
Example: Circumstellar disk with a very small opening angle ^ (see 
Fig. 1, models 2D(a) and 2D(b)) around a pointlike star: Consider- 
ing only those photons being emitted into the disk, the number of 
photons required to achieve a certain accuracy is reduced by a factor 
of 1/ sin^. In contrast to Eq. 2, the angles of emission are now given 

by 

cos = (-1 + 2^6) ■sin(0, (j)E = 2nZj. (26) 
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(2) Taking into account the symmetry of the considered model, the observing 
concept described in Sect. 4 can be modified in order to combine observ- 
ing planes to a larger (efi^ective) solid angle fl. While in the case of ID 
symmetry the angle u can be chosen to be 180°, in models with the z 
axis as an axis of rotation symmetry, all planes within the angular region 
{9o — u) < 9 < [Oq+uj) are combined, whereby is the angle between the 
z axis and the normal vector of the observing plane. The corresponding 
(effective) solid angle VL which is required to normalize the flux/intensity 
then amounts to Q = 27r[cos(6'o — a;) — cos(6'o + a;)]. 

(3) MC3D is prepared to take into account initial dust temperatures Tg > 0, 
such as in active accretion disks around young stellar objects. Therefore, 
the following modifications to the self-consistent RT process have to be 
applied: 

• Iterative reemission: The initial dust radiation field resulting from the 
initial dust temperature and dust density distribution (and dust grain 
optical properties) is added to the absorbed stellar radiation field in 
Eq. 15. 

• Immediate concept: The dust temperatures are set before the stellar 
RT is started. Thus, the absorbed stellar photons provide an additional 
increase of the dust grain temperature. 

(4) MC3D is equipped with a raytracer being optimized for the ID and 2D 
models shown in Fig. 1. It allows derivation of images and the SED 
resulting from the dust reemission radiation. The basing assumption 
that scattering processes are negligible in the corresponding mid-infrared- 
millimeter wavelength range is very well fulfilled for most of the astro- 
physical applications targeted by MC3D. 

Along rays, starting opposite to the observer's location (with respect 
to the model space), the contribution to the monochromatic luminosity 
is integrated taking into account extinction (see Fig. 6[A]). To take into 
account the emission/extinction of the smallest structures but also to 
guarantee a fast integration process, the step width of the numerical 
integration along a ray is adapted to the ESC grid, using a fraction of the 
ESC of the current integration point. The starting points of the rays are 
located in centrosymmetric rings (see Fig. 6[B]), whereby the inner/outer 
radii of these rings correspond to the radial subdivision of the model 
space into ESCs. The area whithin the inner radius of the models (where 
no subdivision of the model into ESCs is provided) is covered by rings 
with the inner/outer radius difference equal to that of the innermost 
ESC(s). Furthermore, the rings are subdivided in the azimuthal direction 
whereby the step width is chosen to achieve a similar resolution as in 
the radial direction. To increase the spatial resolution of the maps, each 
ring can be subdivided in a certain number of rings, and consequently 
smaller azimuthal subdivisions. To project the monochromatic luminosity 
determined in each of the resulting segments (each represented by a single 
ray) as shown in Fig. 6[B] and [C] on the matrix- like observing map, they 
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are subdivided in radial and azimuthal direction again. 
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A Weighted mean optical peirameters of dust grain mixtures 



In case of non-self-consistent radiative transfer (RT) in a dust grain mixture, 
the numerical effort (ultimately, the CPU time and RAM requirement) may be 
substantially decreased by assuming weighted mean values of those parameters 
which describe the interaction of the electromagnetic field with the dust grains. 
The weight Wj (a) which represents the contribution of a certain component of 
the dust grain mixture results from the abundance of the ith material of given 
dust grain number density (assuming tid chemically different dust species) and 
the size distribution of the respective material: 

/ wj{a)da^l. (A.l) 



The Stokes parameters as well as the extinction, absorption, and scattering 
cross section (Cext, C'abs; C'sca) are additive. Therefore, the representative values 
in case of a dust grain mixture can be derived on the basis of their weighted 
contributions (see, e.g., Martin 1978, Sole 1980): 

(Cext) = E / • Cext,- (a) da, (A.2) 

i=j n ■ 

TlD "max 

(Cabs) = J2 J • Cabs,- (a) da, (A.3) 

i=j n ■ 

(Csca) = E y • Csca,-(a) da, (A.4) 



and 



rtD "-ma" 

('^) = X] / '^A^ ' Sj{a)da, (A.5) 



l=J n 



where S is the Miiller matrix which is used to describe the modification of 
the photon's Stokes vector due to the interaction of a photon with the scat- 
tering/absorbing medium (dust grains; see Bickel & Bailey 1985, Bohren & 
Huffman 1983). For the albedo A follows: 

_ x:::;:;r • c., (») • -4,(.) da _ (c^^a) 

^ ^ Sr=°i ^Ao^) ■ ^ext, (a) da (Cext) ■ ^ ■ ^ 
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B Mean free path length 



(1) Decrease of the intensity of a plane wave due to extinction, described by 
Bouguer-Lambert's law 



h 



(B.l) 



(/o ■ ■ .original intensity, 1\ . . . intensity after passing a medium with the 
optical depth r). 
(2) Let Text be the optical depth due to dust extinction: 



(B.2) 



(C'extj is the wavelength-dependent extinction cross section of a single 
dust grain of the species number s is the geometric free path length, 

is the number density of dust particles of that species) 
assumption: Cextj ■ Pj is constant along the geometrical distance s 
(3) Probability /(Text)drext for a photon to pass a distance corresponding to 
the optical depth Text: 



/(Text)dText 



Ip-e '^ext-dText 

/•so : 

J /o-e ■^ext dText 



= e 



-Text 



•dr, 



ext 



(B.3) 



(4) Resulting distribution of the optical depth F(rext): 



Fir. 



ext J 



(B.4) 



(5) Determination of the free optical path langth by substitution of -F(Text) 
by a random number uniformly distributed in the intervall [0, 1]: 



Text = -ln(l - Z) 



(B.5) 



C Scattering matrix 



The scattering Matrix (special Miiller matrix) for the description of the mod- 
ification of the Stokes vector due to the interaction of a weighted photon with 
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a spherical, homogeneous dust particle (Mie scattering) has the form 



^ s^,{e) S^2{0) ^ 

Si2{0) SniO) 

Sssie) Ss4{0) 
\ -Su{e) S2,{e) ) 



(C.l) 



where 



Snie) = li\s^{e)f + \s,{e)f + \Ss{e)f + \s,{e)f) 
s,sie) = Re{s,id)s*id) + s,id)siid)} 

S3,{9) = Re{S2{9)Sm + S^{9)SI{9)} 



51 — J2^=i „(^+i) (O'n'^n + bjiTn), 

52 = E^=l ^gjry («nTn + ftnTTn), (C.3) 
^3 = ^4 = 



■tl]'„(mx)C,n(x)-m^n{rnx)C,!a{x) ' 

(C.4) 

^ _ iTT-i'n (jnx)tpn {x)-ip„ {mx)tp!^ (x) 
" m'il!'„{mx)Cnix)—'4>n(mx)C^{x) 



7rn{cos9) = ^^P^{cos9), 

(C.5) 

T„(cOS^) = |P„l(cOSe) . 

P^(cos6') are associated Legendre functions, ipni^) und Cni^) are Riccati- 
Bessel functions, m is the relative refraction index, and 

2% 

X = — ■ \\mo\\ ■ uk (C.6) 
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is the size parameter. The quantity mo is the refractive index of the surround- 
ing medium|^ A is the wavelength of the infalhng radiation and ok is the dust 
particle radius (see Bohren & Huffman (]T9|)). 

The numerical formulation of the scattering mechanism for spherical, homoge- 
nous particle is based on the algorithm presented by Fischer 



In case of Thomson scattering, the scattering matrix elements are wavelength- 
independent and can be written as 

Sn{e) = S22ie) = {cos\9) + l)/2 
Su{e) = S2ii9) = {cos\9)~l)/2 

= S,,i9) = cosie) (C.7) 

Sl3 = S^i = 5*23 = 5*32 = 
5'l4 = 5*24 = >S'34 = 5*43 = 



(see Bohren & Huffman (|T9D). 

The scattering cross section of an electron in this scattering regime is given 
by 



(e is the elementary charge, Eq is the dielectric constant of free space; see, e.g., 
Musiol et al. (||)). 



D Reemission from ESCs 



In the following the equations for the (random) determination of the reemis- 
sion point within an ESC for the different models shown in Fig. 1 are given. 
According to the symmetry of the particular model, spherical {r,9,(f)), cylin- 
drical {r,(j),z), or cartesian coordinates {x,y,z) are used. 

Definition: Let e be e {x,y,z,r,(j),9}. Then, Si and Sj+i are the inner and outer 
boundary of the ESC along the coordinate e, whereby Sj+i > Si. 



^ Since MC3D was developed to model astrophysical objects for which the sur- 
rounding medium can be assumed to be a perfect vacuum, ||mo||=l. 
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Model: ID 

Definition ranges: 
-Rin < r < i?ou 
0° < ^ < 180° 
0° < (/)< 360° 



cos^ = 1 - 2Z2, 

(f) = 277^3 



(D.l) 
(D.2) 
(D.3) 



Model: 2D (a) 

Definition ranges: 
Rm<r < Ron 

C<9< 180° - ^ 
0° < (^)< 360° 



cos9 = (l-2Z2)sin^ 
(f) = 27rZ3 



(D.4) 
(D.5) 
(D.6) 



Model: 2D(b) 

Definition ranges: 

Rm<r < Ron 

^<e < 180° - c 
0° <d)< 360° 



cos ^ = COS 9i+i + ^2(003 6i — cos ^i+i) 

= 271^3 



(D.7) 
(D.8) 
(D.9) 



Model: 3D (a) 

Definition ranges: 

< ,x < /c 

0<y<lc 

Q< z<L 



X — Xi + Zii^Xi^i — Xi) 

y ^ Vi + Z2{yi+i - Vi) 

z = Zi + Z2{zi^i — Zi) 



(D.IO) 
(D.ll) 

(D.12) 



Model: 3D(b) 

Definition ranges: 
-Rin < r < i?ou 
0° < e < 180° 
0° < (/)< 360° 



- = Z^{ri^^ - r^) + r. 

= 0i + ^2(0i+l - 0i) 

^ = + Z-i{ziJ^\ — Zi) 



(D.13) 
(D.14) 
(D.15) 
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Model 3D(c) 

Definition ranges: 
-Rin < r < i?ou 
0° < ^ < 180° 
0° < (/)< 360° 



cos 6* = cos + Z2{cos9i — cos9i+i] 

(f)^(f)i + Z3{(f)i+i - (f)i) 



(D.16) 
(D.17) 
(D.18) 
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>150K 



100 AU 




Fig. 4. Illustration of the high efficiency of the continuous absorption concept. While 
the temperature distribution shown in the left column was simulated without this 
concept, it was applied to achieve the results shown in the right column. Model: 
Density distribution of the model for the circumstellar disk around the T Tauri 
star HH30 (from Cotera et al. 2001, Wood et al. 2002) but with only 10"^ of ^he 
mass assumed there. The disk is characterized by an optically thick midplane and 
an optically very thin "atmosphere" above the disk. The statistical noise of the 
temperature values is drastically smaller both in the regions of high and low optical 
depth if the continuous absorption concept is applied (the same number of 10^ 
weighted photons was used in both simulations). 
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Fig. 5. [A] Random walk and projection of the weighted photons. [B] Definition 
of the sohd angle fi. The planes 1 and 2 are co-added with the infinite number of 
planes within the solid angle 0,. 




Fig. 6. Raytracer. [a] Projection of the monochromatic luminosity on the observing 
plane, [b] Configuration of the ray location on the observing plane, [c] Projection 
of each mapped region (represented by single ray) onto the matrix-like observing 
plane. The squares 1-4 symbolize single pixels. 
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